坐标转换之高斯投影转大地坐标（C#）
###以西安80坐标举例###
1、大地原点位于陕西省泾阳县永乐镇，椭球参数采用IUG 1975年大会推荐的参数，基本参数为：
2、解算过程包括：
高斯反算和子午线弧长反算
高斯反算需要用到底点纬度，由子午线弧长反算得到，中间涉及到迭代计算；
3、代码
1）Coordin类
1 class Coordin 2 { 3 4 5 //西安80 6 //高斯投影平面坐标值 7 public double x_XA; 8 public double y_XA; 9 public double Hg_XA; 10 11 //大地坐标值(角度°) 12 public double B_XA_Angle; 13 public double L_XA_Angle; 14 15 //大地坐标值(弧长) 16 public double B_XA; 17 public double L_XA; 18 public double H_XA; 19 20 //空间直角坐标值 21 public double X_XA; 22 public double Y_XA; 23 public double Z_XA; 24 }
2）Ellips类
1 //常用椭球参数 2 public class Ellips 3 { 4 //编号 为80时指西安80所用椭球，类推 5 public int num; 6 7 public double a;//长半轴a(m) 8 public double b;//长半轴b(m) 9 public double f;//扁率f 10 public double e1;//第一偏心率e1 11 public double e2;//第二偏心率e2 12 public double LevOrg85;//85水准原点(m) 13 public double L0;//中央子午线(度) 14 15 public class Ellips75_XA80 16 { 17 public const int num = 80; 18 //1975国际椭球_参心 19 public const double a = 6378140.000;//长半轴a(m) 20 public const double b = 6356755.288;//长半轴b(m) 21 public const double f = 1 / 298.257;//扁率f 22 public const double e1 = 0.006694385;//第一偏心率e1 23 public const double e2 = 0.006739502;//第二偏心率e2 24 public const double LevOrg85 = 72.2604;//85水准原点(m) 25 26 } 27 }
3）弧长反算
1 //子午线弧长反解 2 public double LengthInverse(double x_in, double a, double e1) 3 { 4 //基本常量m0,m2,m4,m6,m8 5 double[] mConst = new double[9]; 6 //基本常量a0,a2,a4,a6,a8 7 double[] aConst = new double[9]; 8 9 //计算常量m0,m2,m4,m6,m8 10 mConst[0] = a * (1 - e1 * e1); 11 mConst[2] = 3 * e1 * e1 * mConst[0] / 2.0; 12 mConst[4] = 5 * e1 * e1 * mConst[2] / 4.0; 13 mConst[6] = 7 * e1 * e1 * mConst[4] / 6.0; 14 mConst[8] = 9 * e1 * e1 * mConst[6] / 8.0; 15 //计算常量a0,a2,a4,a6,a8 16 aConst[0] = mConst[0] + mConst[2] / 2 + 3 * mConst[4] / 8
+ 5 * mConst[6] / 16 + 35 * mConst[8] / 128; 17 aConst[2] = mConst[2] / 2 + mConst[4] / 2 + 15 * mConst[6] / 32
+ 7 * mConst[8] / 16; 18 aConst[4] = mConst[4] / 8 + 3 * mConst[6] / 16 + 7 * mConst[8] / 32; 19 aConst[6] = mConst[6] / 32 + mConst[8] / 16; 20 aConst[8] = mConst[8] / 128; 21 22 //设置初始底点纬度B0 23 //进行迭代运算 24 double B0 = x_in / a; 25 bool Judge = true; 26 while (Judge) 27 { 28 double F = -aConst[2] * Math.Sin(2 * B0) / 2
+ aConst[4] * Math.Sin(4 * B0) / 4
- aConst[6] * Math.Sin(6 * B0) / 6
+ aConst[8] * Math.Sin(8 * B0) / 8; 29 double B1 = (x_in - F) / a; 30 if (Math.Abs(B0 - B1) < 1.0E-6) 31 { 32 Judge = false; 33 } 34 else 35 { 36 B0 = B1; 37 } 38 } 39 40 return B0; 41 }
4）高斯反算
1 //高斯反算(xyHg->BLH) 2 public void GaussBacCalculation(Coordin coordin, Ellips ellips) 3 { 4 //初始化椭球参数 5 int num = ellips.num; 6 double a = ellips.a; 7 double e1 = ellips.e1; 8 double e2 = ellips.e2; 9 double L0 = ellips.L0; 10 double LevOrg85 = ellips.LevOrg85; 11 //初始化为0 12 double x = 0; 13 double y = 0; 14 double Hg = 0; 15 //根据椭球判断 16 switch(num) 17 { 18 case 54: 19 //北京54高斯投影平面坐标 20 21 break; 22 23 case 80: 24 //西安80高斯投影平面坐标 25 x = coordin.x_XA; 26 y = coordin.y_XA; 27 Hg = coordin.Hg_XA; 28 break; 29 30 } 31 32 33 34 //(弧长)解算底点纬度Bf 35 double Bf = LengthInverse(x, a, e1); 36 //中间变量Wf 37 double Wf = Math.Sqrt(1 - Math.Pow(e1, 2) * Math.Sin(Bf)); 38 //子午圈曲率半径Nf 39 double Nf = a / Wf; 40 //其余变量 41 double Mf = a * (1 - Math.Pow(e1, 2)) / (Math.Pow(Wf, 3)); 42 double tf = Math.Tan(Bf); 43 double pf = e2 * Math.Cos(Bf); 44 //修正横坐标y值 45 y -= 500000; 46 47 //(弧长)经度L 48 double l = y / (Nf * Math.Cos(Bf)) 49 - Math.Pow(y, 3) * (
1 + 2 * Math.Pow(tf, 2) + Math.Pow(pf, 2)
) / (6 * Math.Pow(Nf, 3) * Math.Cos(Bf)) 50 + Math.Pow(y, 5) * (
5 + 28 * Math.Pow(tf, 2) + 24 * Math.Pow(tf, 4)
+ 6 * Math.Pow(pf, 2) + 8 * Math.Pow(pf, 2) * Math.Pow(tf, 2)
) / (120 * Math.Pow(Nf, 5) * Math.Cos(Bf)); 51 double L = l + L0 * rPI; 52 //(弧长)维度B 53 double B = Bf - tf * Math.Pow(y, 2) / (2 * Mf * Nf)
+ tf * Math.Pow(y, 4) * (
5 + 3 * Math.Pow(tf, 2) + Math.Pow(pf, 2)
- 9 * Math.Pow(pf, 2) * Math.Pow(tf, 2)
) / (24 * Mf * Math.Pow(Nf, 3)) 54 - tf * Math.Pow(y, 6) * (
61 + 90 * Math.Pow(tf, 2) + 45 * Math.Pow(tf, 4)
) / (720 * Mf * Math.Pow(Nf, 5)); 55 //(不考虑高程异常)正常高转大地高 56 double H = Hg + LevOrg85; 57 58 59 switch (num) 60 { 61 case 54: 62 //北京54大地坐标 63 64 break; 65 66 case 80: 67 //西安80大地坐标 68 coordin.L_XA = L; 69 coordin.B_XA = B; 70 coordin.H_XA = H; 71 break; 72 } 73 74 75 }
参考博客及文档
