批量计算有利于提高工作/学习效率，本文以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();}}}
