C#实现批量高斯投影正算、反算

​欢迎关注博主的微信公众号：跟着MO学习GIS

批量计算有利于提高工作/学习效率，本文以EPSG提供《Coordinate Conversions and Transformations including Formulas》
的高斯投影正算、反算算法写成c#代码为例。。

// 高斯投影反算，将高斯坐标反算出经纬度坐标

class Gausstojw

    {

        private double atanh(double z)

        {

            return 0.5 * Math.Log((1 + z) / (1 - z));

        }


        private double asinh(double z)

        {

            return Math.Log(z + Math.Sqrt(z * z + 1));

        }



        private double ψ = -9999;

        private double lambda = -9999;


        public double Longitude

        {

            get

            {

                if (lambda == -9999)

                {


                    GuassToGeo(_a, _b, 1 / _invf, _FE, _FN, _X, _Y, _CM, _k);

                }

                return lambda;

            }

        }


        private double _a = 6378245.0;

        private double _b = (298.300000000000010000 - 1) / 298.300000000000010000 * 6378245.0;

        private double _invf = 298.300000000000010000;

        private double _FE = 500000;

        private double _FN = 0;

        private double _CM = 111;

        private double _k = 1;


        public int Zone

        {

            set

            {

                _CM = value * 6 - 3;

            }

        }


        public double SemiMajor

        {

            set

            {

                _a = value;

            }

            get

            {

                return _a;

            }

        }


        public double Lantitude

        {

            get

            {

                if (lambda == -9999)

                {

                    GuassToGeo(_a, _b, 1 / _invf, _FE, _FN, _X, _Y, _CM, _k);

                }

                return ψ;

            }

        }


        private double _X;

        private double _Y;

        public double X

        {

            set

            {

                _X = value;

            }

        }


        public double Y

        {


            set

            {

                _Y = value;

            }

        }


        private void GuassToGeo(double a, double b, double f, double FE, double FN, double E, double N, double λ0, double k0)

        {

            //double f = (a - b) / a;

            double n = f / (2 - f);

            double e = Math.Sqrt((a * a - b * b)) / a;

            double B = (a / (1 + n)) * (1 + n * n / 4 + n * n * n * n / 64);


            double h1 = (n / 2) - (2.0 / 3) * n * n + (37.0 / 96) * n * n * n - (1.0 / 360) * n * n * n * n;

            double h2 = (1.0 / 48) * n * n + (1.0 / 15) * n * n * n - (437.0 / 1440) * n * n * n * n;

            double h3 = (17.0 / 480) * n * n * n - (37.0 / 840) * n * n * n * n;

            double h4 = (4397.0 / 161280) * n * n * n * n;


            double η = (E - FE) / (B * k0);

            double ζ = (N - FN) / (B * k0);

            double ζ1 = h1 * Math.Sin(2 * ζ) * Math.Cosh(2 * η);

            double ζ2 = h2 * Math.Sin(4 * ζ) * Math.Cosh(4 * η);

            double ζ3 = h3 * Math.Sin(6 * ζ) * Math.Cosh(6 * η);

            double ζ4 = h4 * Math.Sin(8 * ζ) * Math.Cosh(8 * η);

            double ζ0 = ζ - (ζ1 + ζ2 + ζ3 + ζ4);


            double η1 = h1 * Math.Cos(2 * ζ) * Math.Sinh(2 * η);

            double η2 = h2 * Math.Cos(4 * ζ) * Math.Sinh(4 * η);

            double η3 = h3 * Math.Cos(6 * ζ) * Math.Sinh(6 * η);

            double η4 = h4 * Math.Cos(8 * ζ) * Math.Sinh(8 * η);

            double η0 = η - (η1 + η2 + η3 + η4);


            double Β = Math.Asin(Math.Sin(ζ0) / Math.Cosh(η0));

            double Q = asinh(Math.Tan(Β));

            double Q1 = Q + (e * atanh(e * Math.Tanh(Q)));

            for (int i = 1; i < 400; i++)

            {

                Q1 = Q + (e * atanh(e * Math.Tanh(Q1)));

            }


            ψ = Math.Atan(Math.Sinh(Q1)) * 180 / 3.1415926;

            lambda = λ0 + Math.Asin(Math.Tanh(η0) / Math.Cos(Β)) * 180 / 3.1415926;

            //System.Console.WriteLine("wd:{0},jd:{1}", ψ, λ);



        }


//球面坐标经纬度投影成高斯坐标

class JWtoGauss

    {

        private double E = -9999;

        private double N = -9999;

        private double a3 = -9999;

        public double _a3

        {

            set

            {

                if (a3 == -9999)

                {

                    GEOToGuass(_a, _b, _FE, _FN, _E, _N, _CM, _Nψ0, k);

                }

            }

            get { return a3; }

        }


        private double a2 = -9999;

        public double _a2

        {

            set

            {

                if (a2 == -9999)

                {

                    GEOToGuass(_a, _b, _FE, _FN, _E, _N, _CM, _Nψ0, k);

                }

            }

            get { return a2; }

        }

        private double _a = 6378245.0;

        public double longR

        {

            set { _a = value; }

            get { return _a; }

        }

        private double _b = 6356863.01877304730;

        public double shortR

        {

            set { _b = value; }

            get { return _b; }

        }

        private double _FE = 500000.0;

        private double _FN = 0.0;



        private double _CM = 111.0 / (180 / 3.1415926);//注意顺序括号在后！

        public double jisuancenter

        {

            set { _CM = CenterLonNum6(value); }


        }

        private double _Nψ0 = 0.0;

        private double k = 1.0;

        public double X

        {

            get

            {

                if (E == -9999)

                {

                    GEOToGuass(_a, _b, _FE, _FN, _E, _N, _CM, _Nψ0, k);

                }

                return E;

            }

        }

        public double Y

        {

            get

            {

                if (N == -9999)

                {

                    GEOToGuass(_a, _b, _FE, _FN, _E, _N, _CM, _Nψ0, k);

                }

                return N;

            }

        }

        private double _E;

        public double _EX

        {

            set { _E = value / 57.29578; }

        }

        private double _N;

        public double _NY

        {

            set { _N = value / 57.29578; }

        }

        private double atanh(double z)

        {

            return 0.5 * Math.Log((1 + z) / (1 - z));

        }

        private double asinh(double z)

        {

            return Math.Log(z + Math.Sqrt(z * z + 1));

        }

        private double CenterLonNum6(double Eλ)

        {

            int a = (int)((Eλ * 180 / 3.1415926) / 6);

            double center = ((a + 1) * 6 - 3);

            return center / (180 / 3.1415926);


        }

        private double CenterLonNum3(double Eλ)

        {

            //if (Eλ < 1.5)

            //{

            //    //jc=jd;

            //    center = Eλ;

            //}

            double n = Eλ - 1.5;

            int center = ((int)(n / 3) + 1) * 3;

            //int center=n*3;


            return center * 1.0;

        }

        private void GEOToGuass(double a, double b, double FE, double FN, double Eλ, double Nψ, double Eλ0, double Nψ0, double k0)

        {

            double f = (a - b) / a;

            double e = Math.Sqrt((a * a - b * b)) / a;

            double n = f / (2 - f);


            double B = (a / (1 + n)) * (1 + n * n / 4 + n * n * n * n / 64);

            double h1 = n / 2 - ((2.0 / 3) * n * n) + ((5.0 / 16) * n * n * n) + ((41.0 / 180) * n * n * n * n);

            double h2 = ((13.0 / 48) * n * n) - ((3.0 / 5) * n * n * n) + ((557.0 / 1440) * n * n * n * n);

            double h3 = ((61.0 / 240) * n * n * n) - ((103.0 / 140) * n * n * n * n);

            double h4 = (49561.0 / 161280) * n * n * n * n;

            //double ψ=0.0;


            //double Q0 = asinh(Math.Tan(ψ))-(e*atanh(e*Math.Sin(ψ)));

            //double β0 = atanh(Math.Sinh(Q0));

            //double ξq0 = Math.Asin(Math.Sin(β0));

            double Q = asinh(Math.Tan(Nψ)) - (e * atanh(e * Math.Sin(Nψ)));

            double β = Math.Atan(Math.Sinh(Q));


            double η0 = atanh(Math.Cos(β) * Math.Sin(Eλ - Eλ0));

            //double η = (E - FE) / (B * k0);

            //double ζ = (N - FN) / (B * k0);

            double ζ0 = Math.Asin(Math.Sin(β) * Math.Cosh(η0));


            double ζ1 = h1 * Math.Sin(2 * ζ0) * Math.Cosh(2 * η0);

            double ζ2 = h2 * Math.Sin(4 * ζ0) * Math.Cosh(4 * η0);

            double ζ3 = h3 * Math.Sin(6 * ζ0) * Math.Cosh(6 * η0);

            double ζ4 = h4 * Math.Sin(8 * ζ0) * Math.Cosh(8 * η0);

            double ζ = ζ0 + ζ1 + ζ2 + ζ3 + ζ4;


            double η1 = h1 * Math.Cos(2 * ζ0) * Math.Sinh(2 * η0);

            double η2 = h2 * Math.Cos(4 * ζ0) * Math.Sinh(4 * η0);

            double η3 = h3 * Math.Cos(6 * ζ0) * Math.Sinh(6 * η0);

            double η4 = h4 * Math.Cos(8 * ζ0) * Math.Sinh(8 * η0);

            double η = η0 + η1 + η2 + η3 + η4;


            E = FE + k0 * B * η;

            N = FN + k0 * B * ζ;

            a3 = 4495668.7826 - N;

            a2 = 542576.2338 - E;



        }

    }


//主程序

class Program

    {

        static void Main(string[] args)

        {

            

            //批量投影算法

            //打开两个文本文件，一个读取数据，一个用来存放数据

            string file = @"C:\Users\suns\Desktop\POI1.txt";

            string W = @"C:\Users\suns\Desktop\1.txt";

            FileStream f = File.Open(W, FileMode.Open, FileAccess.ReadWrite);//用于存放计算后的数据以可读可写的权限打开

            StreamReader r = File.OpenText(file);//StreamReader 以打开

            StreamWriter w = new StreamWriter(f);//向1.txt写入东西的声明

            string str = "";

            int num = 0;

           //一个循环读取文本中所有存放数据的行，调用上面的类进行批量计算然后写入另一个文本中

            while ((str = r.ReadLine()) != null)

            {

                JWtoGauss a1 = new JWtoGauss();

                num++;

                string[] fs = str.Split(new char[] { ',' });//行的数据格式:xxx,yyy。以“，”将它们分成两个部分，用一个字符串数组存起来

                //都将x,y转成双精度型数字

                Double a = Convert.ToDouble(fs[0]);

                Double b = Convert.ToDouble(fs[1]);

                a1._EX = a;

                a1._NY = b;

                w.WriteLine(a1.X + "," + a1.Y);

            }

            w.Flush();

            w.Close();

            f.Close();

            Console.WriteLine(num);

            Console.Read();

————————————————————————————————————————————————————————————————————


           //批量高斯投影反算算法

            static void Main(string[] args)

            {

                

                string file = @"C:\Users\suns\Desktop\1.txt";

                string W = @"C:\Users\suns\Desktop\2.txt";

                FileStream f = File.Open(W, FileMode.Open, FileAccess.ReadWrite);

                StreamReader r = File.OpenText(file);

                StreamWriter w = new StreamWriter(f);

                string str = "";

                int num = 0;

                while ((str = r.ReadLine()) != null)

                {

                    Gausstojw a1 = new Gausstojw();

                   // jwToGuass a1 = new jwToGuass();

                    num++;

                    //char i = ','; 

                    string[] fs = str.Split(new char[] { ',' });

                    Double a = Convert.ToDouble(fs[0]);

                    Double b = Convert.ToDouble(fs[1]);

                    a1.X = a;

                    a1.Y= b;

                    a1.Zone=20;

                    w.WriteLine(a1.Longitude + "," + a1.Lantitude);

                }

                w.Flush();

                w.Close();

                f.Close();

                Console.WriteLine(num);

                Console.Read();


            }



        }


    }