2坐标格式转换 ∈ C# 编程笔记

【日志】
2020/6/25 开了这篇博客，功能还算满意

一、原理 - 并不在这里

原理部分请看博客1：

https://blog.csdn.net/Gou_Hailong/article/details/106938497

二、源码

源码结构：

	class zbgsTransfor
    {
   
        /*辅助函数 - 用于求d的n次方,因为用Math.Pow()太累赘
         */ 
        public static double P(double d, double n)
        {
   
            return Math.Pow(d, n);
        }
        /*辅助函数 - 用于计算椭球参数
         * 输入：a 长轴, b 短轴,
         * 输出：tq 椭球参数数组
         *        0 a, 1 b, 2 c, 3 e
         *        4 e1, 5 alpha
         */
        static double[] Get_TQ(double a, double b)
        {
   
            double[] tq = new double[6];
            double c = a * a / b,
                e = Math.Sqrt(a * a - b * b) / a,
                e1 = Math.Sqrt(a * a - b * b) / b,
                alpha = (a - b) / a;
            tq[0] = a; tq[1] = b; tq[2] = c;
            tq[3] = e; tq[4] = e1; tq[5] = alpha;
            return tq;
        }
        /*此函数用来将xyz 转换成blh
         * 输入：pt 待转换的点数组, mod 选择哪个椭球 1 75 2 克
         * 输出：转换成功的点数组
         * 注：默认WGS84椭球参数，可改
         */
        public static Point[] Xyz2Blh(Point[] pt, int mod=1)
        {
   
            int n = pt.Length;
            Point[] Pt = new Point[n];
            double a = 0, b = 0, e2, e1;
            if (mod == 1)
            {
   //1975年椭球
                a = 6378140;
                b = 6356755.2881575287;
            }
            else if (mod == 2)
            {
   //克拉索夫斯基椭球
                a = 6378245.0;
                b = 6356863.0187730473;
            }
            double[] TQ = Get_TQ(a, b);
            e1 = TQ[3];
            e2 = TQ[4];
            for(int i=0;i<n;i++)
            {
   
                Point p=pt[i];            
                double B, B0 = 0, L, H, N;
                L = Math.Atan(p.y / p.x);
                B = Math.Atan(p.z / Math.Sqrt(P(p.x, 2) + P(p.y, 2)));
                while (Math.Abs(B0 - B) > 1e-9)
                {
   
                    B0 = B;
                    N = a / Math.Sqrt(1 - P(e1 * Math.Sin(B), 2));
                    B = Math.Atan((p.z + N * P(e1, 2) * Math.Sin(B)) / Math.Sqrt(P(p.x, 2) + P(p.y, 2)));
                }
                N = a / Math.Sqrt(1 - P(e1 * Math.Sin(B), 2));
                H = Math.Sqrt(P(p.x, 2) + P(p.y, 2)) / Math.Cos(B) - N;
                Pt[i] = new Point(B, L, H, p.name);
            }
            return Pt;
        }
        /*此函数用来将blh 转换成xyz
         * bl都是弧度！！
         * 输入：待转换的点数组
         * 输出：转换成功的点数组
         * 注：默认WGS84椭球参数，可改
         */
        public static Point[] Blh2Xyz(Point[] pt, int mod = 1)
        {
   
            int n = pt.Length;
            Point[] Pt = new Point[n];
            double a = 0, b = 0, e2, e1;
            if (mod == 1)
            {
   //1975年椭球
                a = 6378140;
                b = 6356755.2881575287;
            }
            else if (mod == 2)
            {
   //克拉索夫斯基椭球
                a = 6378245.0;
                b = 6356863.0187730473;
            }
            double[] TQ = Get_TQ(a, b);
            e1 = TQ[3];
            e2 = TQ[4];
            for(int i=0;i<n;i++)
            {
   
                Point p=pt[i];           
                double N = a / Math.Sqrt(1 - Math.Pow(e1 * Math.Sin(p.x), 2));
                double x, y, z;
                x = (N + p.z) * Math.Cos(p.x) * Math.Cos(p.y);
                y = (N + p.z) * Math.Cos(p.x) * Math.Sin(p.y)